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Abstract. We derive a systematic, multiple time-scale perturbation expansion for 
the work distribution in isothermal quasi-static Langevin processes. To first order we 
find a Gaussian distribution reproducing the result of Speck and Seifert [Phys. Rev. E 
70, 066112 (2004)]. Scrutinizing the applicability of perturbation theory we then show 
that, irrespective of time-scale separation, the expansion breaks down when applied to 
untypical work values from the tails of the distribution. We thus reconcile the result 
of Speck and Seifert with apparently conflicting exact expressions for the asymptotics 
of work distributions in special systems and with an intuitive argument building on 
the central limit theorem. 
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1. Introduction 

Thermodynamics of systems so small that typical changes of their energies are of the 
order of ksT is concerned with distributions of thermodynamic quantities rather than 
with their averages [U Ej [3] . These distributions fulfill a number of exact and general 
relations which are now commonly referred to as fluctuation theorems. At least two 
features of these theorems are truly remarkable: First, they hold for systems driven 
(almost) arbitrarily far from equilibrium, and second, they are very sensitive to the tails 
of the respective probability distributions 

The concept of work is a cornerstone of thermodynamics. Accordingly, the 
distribution of work is of central importance in stochastic thermodynamics. If the 
parameters of a system coupled to a heat bath are changed very slowly one expects 
that the system stays in equilibrium during the whole process. In this case the work 
performed or consumed in the transition is not fluctuating and equals the difference in 
free energy between final and initial state of the system. If, on the other hand, the 
system is driven violently it will pass through intermediate non-equilibrium states. In 
this case a fluctuating fraction of dissipative work adds to the free-energy difference. The 
resulting distribution of work is, in general, non-universal displaying features specific 
for the system and the process at hand. 

In an attempt to establish some universality even for non-equilibrium processes 
Speck and Seifert investigated the work distribution for isothermal stochastic processes 
with small but non-zero driving [5]. By using a projection-operator technique to derive 
approximate solutions of the underlying Fokker-Planck equation [6] they found that in 
these cases the work distribution must be Gaussian to leading order. A simple and 
intuitive argument put forward in [7] relates this result to the central limit theorem. 
In a quasi-static process the system has time to relax to its instantaneous equilibrium 
distribution between changes of the external parameters. Therefore, the total work 
becomes a sum of many independent contributions and by virtue of the central limit 
theorem the resulting distribution has to converge to a Gaussian. 

On the other hand, the determination of the exact asymptotics of the work 
distribution for a simple model system found an exponential tail for any driving [8]. 
Hence, the Gaussian character of the work distribution of quasi-static processes does 
not extend to the rare realizations. This is again in accordance with the central limit 
theorem which prescribes a Gaussian for the central part of the distribution leaving the 
tails unspecified [9]. From the systematic analysis in [5] it appeared, however, that time 
scale separation alone is sufficient to derive a Gaussian form of the work distribution. 

In the present note we analyze the work distribution in quasi-static stochastic 
processes by employing a multiple-scale perturbation expansion. The expansion 
parameter quantifies the time-scale separation in the problem. We first rederive the 
results of Speck and Seifert and then show that, irrespective of time-scale separation, 
the expansion fails for the tails of the work distribution. A similar argument applies to 
the projection operator technique used by Speck and Seifert. We illustrate the point 
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with numerical simulations for a model system and also indicate how the expansion may 
be extended to higher orders. 

2. The model 

As a special yet representative example we consider the overdamped dynamics of a 
degree of freedom x in a time-dependent potential V{x^ A). The protocol A(t) specifies 
the time-dependence of the potential and evolves from A(0) = at the beginning of 
the process to A(t/) = 1 at its end. Upon redefinition of A we may always achieve a 
linear protocol, A = t/tf. For notational simplicity we restrict ourselves to scalar x, 
generalization to higher dimensions being straightforward. The dynamics is given by 

9tx(t) = -/i9,.r(x,A) + C(t), (1) 
where denotes the mobility and is a standard white-noise source with correlation 

mm) = '^-jKt-t'). (2) 

Here (5 denotes the inverse temperature of the heat bath. The work W performed along 
a trajectory x{t) is given by [TUl E] 

W[x{-)]= rdtXdxV{x{t),X{t)) , (3) 

hence 

dtW = XdxVix,X), (4) 

where the dot denotes the time derivative. 

The joint probability distribution p{x, W, t) describing the coupled stochastic 
evolution of x and W as given by ([1]) and (jlj) therefore fulfills the Fokker-Planck equation 

dtp = fid, [{d.V)p^ - X dxV dwP +^dlp. (5) 

It is convenient to write this equation in the form 

dtp = (L, + Lf ) p (6) 



with 



and 



fid4d,V{x,X) + ^d,) (7) 



if ■.= -Xd,Vix,X)dw . (8) 

For each fixed value A of the protocol parameter Lx has a right eigenvector with 
eigenvalue zero given by the equilibrium distribution 

f,(x) := e'^(^^-^(-'^)) (9) 
with the free energy 

f.:^-iln/..e-<-). (10) 
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The corresponding left eigenvector is given by f{x) = 1. We assume that the process 
starts in equilibrium. According to (|3]) the initial condition for (|5]) therefore reads 

v{x,W,Q) = h{x)5{W) . (11) 

Our aim is to characterize the marginal work distribution 

P{W,t) := j dx p{x,W,t) (12) 

in the quasi-static limit, i.e. when A changes sufficiently slowly. 

3. Multiple time-scale perturbation theory 

The concept of a quasi-static process rests on time-scale separation. For fixed A the 
relaxation of the marginal distribution f{x,t) := J dWp{x,W,t) to its equilibrium 
form, f\{x), happens with a characteristic time Trei, which, in general, depends on A. 
We call a process quasi-static, if 

A(t)rrci<l VtG(0,t/). (13) 

Physically, this means that X(t) hardly changes over times of order Tj-^i or, equivalently, 
that on the time scale on which A evolves f{x, t) is always very near to f\(t){x)- Choosing 
Tj-ei as time unit it is therefore natural to use e := X = 1/tf a.s a. small parameter and to 
investigate the problem perturbatively. 

For e = no work at all is performed and P{W,t) = 5{W) for all t. The form of 
P(W, t) changes qualitatively when e becomes non-zero, i.e. the perturbation is singular 
and a multiple-scale perturbation scheme is called for [HI [12]. Accordingly, we introduce 
an additional slow time variable T := et = X and perform the substitutions 

p{x, W, t) ^ p{x, W, t. A), dt^dt + edx . (14) 

Moreover, since our focus is on the work distribution evolving on the slow time scale 
A, we are not interested in the fast dynamics and neglect the dependence of p{x, W, t, A) 
on t altogether. The Fokker-Planck equation (|6]) then acquires the form 

L^pix, W,X)=e {dx + dxV{x, A) dw) p{x, W, A) (15) 

appropriate for a perturbation expansion. The singular character of the perturbation is 
apparent from the fact that the time derivative is among the small terms. 
We solve ( !T5l) iteratively using the ansatz 

p{x, W, A) = po{x, W, A) + epi{x, W, A) + e^p^ix, W,X) + ... . (16) 

Omitting the t-dependence of p may be interpreted as part of this ansatz. 

Plugging f lT6|) into f lTSj) and matching powers of e we generate a hierarchy of 
equations for the unknown functions p„(x, W, A). To zeroth order we find 

LxPo = (17) 
with the solution 

Po{x,W,X) = fx{x)go{W,X), (18) 
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where go{W, A) is a so far undetermined function constrained only by the normahzation 
condition 



VA . 



(19) 



dWgoiW, A) = 1 
To order e we get 

Lx Pi = {dx + dxV dw) Po = (dxfx) 9o + fx{dx9o + dxV dw9o) ■ (20) 

In order for the perturbation expansion to be well-defined and secular terms being absent 
the r.h.s. of this equation must be orthogonal to the null space of the adjoint operator 

of Lx , i.e. orthogonal to f{x) = 1. We hence impose the solvability condition 



21 







dx 



{dxfx)9o + fx{dx9o + dxV dw9o) 



Now 



(22) 



(23) 



dxdxfxix) = 9x j dxfx{x) = 
and (l22l) translates into 

dx9o = -Aodw9o (24) 

with 

Ao(A) := {dxV)x := J dx fx{x) dxV{x, A) = dxFx . (25) 
Using ( l24l) as well as 

dxfx = /3(Ao - dxV)fx (26) 
as follows from (Q in the r.h.s. of (!20|) we get 

Lxpiix, W, A) = /A(a;)ao(a;, A) [f3go{W, A) - dw9o{W, A)] , (27) 

where 

ao{x,\):=Ao{\)-dxVix,\). (28) 
Eq. ( 1271) may now be solved. Observing 

L-\fx{x)h{x)) = fx{x){Ltr'h{x) (29) 
for any function h{x) [6] we may write the solution in the form 

Pi{x, W, A) = fx{x)bo{x, A) [f3go{W, A) - dw9o{W, A)] + fx{x)9i{W, A)(30) 

where 

bo:={L+y'ao (31) 

denotes the normalizable solution of the equation bo = ao. The second term on 
the r.h.s. of (1301) is the general solution of the homogeneous equation LxPi = 0. By 
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a redefinition of gi we may include any contribution proportional f\{x) from the first 
term into tliis second one. Witliout loss of generality we may therefore require 



dxf{x)h{x)ho{x, A) = (6o)a = . (32) 

At order we have 

L\P2 = (dx + dxVdw)pi ■ (33) 
Using ( 130|) . (123|) . and ( 132|) the solvability condition for this equation takes the form 

dxgi = -AodwQi - Ai{pdwgo - d^go) , (34) 

with 

Ai{X):={dxVbo)x. (35) 

These results are sufficient to derive an equation for P{W, A) correct up to order e. To 
this order it is not necessary to actually solve ( l33l) . 

4. The Gaussian approximation 

From (in]), (US]), dan]), and we find 

P{W, A) = go{W, A) + eg,{W, A) + 0{e^) . (36) 

Therefore 

dxP = dxgo + edxgi + 0{e^) (37) 
which using ( 124)) and ( 134)) gives rise to 

dxP = - Aodwgo - eAodwgi - eAi{l3dwgo - d^go) + Oie"^) 
= - Aodwigo + egi) - eAi{(3dwgo - d^go) + 0{e^) 
= - AodwP - sA.i^dwP - d^^P) + 0{e^) 

= -{Ao + ef3A^)dwP + eA^d^P + 0{e^) . (38) 
The solution of this equation satisfying the initial condition ( ITT)) is a Gaussian 

P(W, A) = ^= exp (-"^-.'^it'n (39) 



with 



and 



a^(A) = 2£ t dX'A^iX') (40) 

^0 



iy(A) = d\' [Ao(A') + e(3A^iX')] = AF(A) + ^a^(A) . (41) 







Eqs. (I59)) - (1H)) reproduce the results of Speck and Seifert [S]. 
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5. Self-consistency of the perturbation expansion 

We next address the question under which conditions our perturbative treatment is self- 
consistent. To this end it is instructive to consider the distribution in x conditioned to 
a given value of W . From ( fT6|) . ( !36!) . as well as from (fTSl) and ( !30|) we find 

p(x,iy,A) 



p{x,\\W) :-- 



P{W,X) 
J\[x) 



= h{x)^ + eho{l3-dw\nP) + 0{e^)\ . (42) 

For perturbation theory to be applicable we hence need 

e|6o(a;,A)(/3-9vKlnP(iy,A))| < 1 . (43) 

Consequently, e <^ 1, as ensured by time-scale separation alone, is not sufficient to make 
the perturbation expansion meaningful. In the present context we are not concerned 
with the singular character of the small-noise limit, /3 — )■ oo, and assume /3 = 0{1). The 
crucial criterion for the self-consistency of our perturbation expansion is therefore 

e|6o<9vFlnP| < 1 . (44) 

From ([39]) we find \dwlnP{W)\ = \W - W\/a^. Hence ([39]) gives a rehable 
approximation of the true work distribution for typical values of W only. For values of 
W sufficiently different from its average (as measured in units of the standard deviation 
aw) P{W, A) may signiffcantly differ from the Gaussian form. In fact, since bo does not 
depend on W there is for every value of £ a threshold Wc such that fl39l) substantially 
deviates from the true work distribution when |iy| ^ Wc- 

From fl42p we realize that the failure of perturbation theory is related to the fact that 
p{x, X\W) may be rather different from f\{x). In fact this makes perfect sense. Typical 
values of W are produced by typical trajectories x(t) the endpoints x(tf) of which provide 
a fair sampling of f\{x). Conditioning the distribution of x to unusual work values, on 
the other hand, implies a strong bias on the weights of the trajectories which generically 
result in substantial deviations of p{x, X\W) from the equilibrium distribution f\{x). To 
illustrate this fact we show in Fig. [1] results obtained by simulating ([1]) for the "breathing 
parabola" 

Vix,X) = ^x' (45) 

with /i = /3 = 1 and tf = 10. The left figure displays the work distribution on a 
logarithmic scale and shows that for most VT-values the distribution is indeed not 
Gaussian. In the right figure we compare the equilibrium distribution fx{x) with 
p{x, X\W) for different intervals of H^- values. As expected, choosing W from the center 
of P{W,X) gives rise to a p{x,X\W) rather similar to fx{x). Very small values of W, 
however, require trajectories x{t) with sufficiently large excursions from x = such 
that p{x, X\W) broadens and finally develops a bimodal shape that is even qualitatively 
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Figure 1. Left: Results from a numerical simulation of Eq. ([T]) for the "breathing 
parabola" (P5|) with /i = /3 = 1 and = 10 together with the perturbative results 
of first and second order. Note the strong deviation between the distribution P{W) 
as obtained from the simulations (blue line) and the Gaussian approximation (j39p 
(red line). Right: Comparison of the final equilibrium distribution /i(a;) with the 
conditional distribution p{x\W) for several intervals of W . All histograms are compiled 
from the simulation of 6 • 10^ trajectories, the conditional distributions are smoothed 
by a spline-interpolation. 



different from the Gaussian equilibrium distribution. Note that for the potential chosen 
Trei = l/(/x(2 — A)) < 1 for the entire process. Together with A = l/tj = 0.1 time-scale 
separation as defined by ( fT3l) is hence well satisfied. 



6. Higher orders 

It is straightforward to extend our perturbation expansion to higher orders. Enforcing 
we rewrite (1331) in the form 



Lx P2 = fx 



aoi/^di - dwgi) + (/3ao6o + dxbo){/3go - dwQo) 
- (Ao6o + ai){l3dwgo - d^go) , 



with 



ai(x,A) ■.= Ai{X)-dxV{x,X)bo{x,X) 
The solution of (l46l) is given by 

P2 = fx 



(46) 
(47) 



boiPgi - dwgi) + bi{/3go - dwgo) - ci{(3dwgo - d^go) + g2 , (48) 
where we have introduced the normalizable solutions 

6i(x, A) := (L+ r\f3aobo + dxbo) and (49) 
Ci(x,A):=(L+)-i(Ao6o + ai). (50) 
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as well as another solution f\{x)g2{W,\) of the homogeneous equation. Similarly to 
(132|) we may assume that 

= {bi)x{l3go - dxgo) - {ci)\{(3dwgo - <9^fl'o) • (51) 
Using (148!) we may explicitly calculate the r.h.s. of the order equation 

LxP3={dx-dxVdw)P2. (52) 
The solvability condition at this order takes the form 

dxg2 = - Aodwg2 - Ai{l3dwgi - d^gi) - A2il3dwgo - d^go) 

+ C2(/39^^?o-9^^?o), (53) 

with 

^(A) := {dxVbi)x and C2(A) := {dxVci)x ■ (54) 

Proceeding in the same way as we did in deriving (l38l) we end up with 

dxP = - (Ao + e(3A, + e^(3A2)dwP + {eA, + e^A^ + /3C2))9^P 

- e^C^d^^P + 0{e^) . (55) 

The last term describes deviations from a Gaussian P{W, A) which hence first show up 
at order e"^. 

As a simple example we have carried through this program for the potential f HSj) . 
The results are 

11 

"2/3(2 -A)' ~ 4/./32(2 - A)3 ' ^^^^ 

and 

3 1 

" ~8/i2/32(2- A)5' " "4/^2/33(2 _A)5 ' ^^^^ 

The left plot in Fig. [1] compares the first and second order approximation of P{W) with 
the histogram resulting from numerical simulations. Using the expressions for Aq and 
Ai in (HOj) and (HTj) gives rise to the Gaussian approximation shown by the red line. The 
data for the second order approximation shown in green were generated by numerically 
solving fl55|) with the (5-function in the initial condition f|TT]) replaced by a very narrow 
Gaussian. One clearly sees that the first order Gaussian approximation fails for most 
of the VT-values. The second order result improves on the first order but is nevertheless 
still far from a satisfactory approximation for most of the distribution. 

Let us also note, that although our perturbation expansion is complementary to 
the projection operator analysis of Speck and Seifert both approaches are, of course, 
related. In particular, conditions ( l32l) and ( ISTl) single the part proportional to fx{x) 
out of the solutions for pi{x, W, A) and P2{x, W, A) respectively. This is analogous to the 
effect of the projection operator IIa introduced by Speck and Seifert [5]. 
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7. The "simplest" example 

A somewhat curious case is given by the so-called "shifted parabola" 

V{x,\) = ^{x-\f (58) 

that is often considered to be the simplest example of a driven Langevin system. In fact 
in this case the distribution of work is Gaussian for any driving speed \ = e [131 E 



P(iy, A) = -jL= exp ('_2_JQ! ) (59) 



with 

W{X) = -\ + ^{e-^^^'' -I) and aU\) = lw{X). (60) 

The perturbative treatment of this case is, however, somewhat subtle. First, we get 
from (125]) Aq = consistent with the fact that the free energy does not change with A. 
Accordingly, qq does not depend on A either (cf. (^^)- The same holds true for 
Ai = l/(/i/3). Consequently, the complete r.h.s. of (134|) is independent of A which 
implies gi ~ A, i.e. generates yet another secular term in the perturbation expansion. 
A way to get rid of this extra complication is to introduce the "super-slow" time scale, 
T, vi&t t + eX + e^T Instead of (El]) we then find 

dxgi = -Aodwgi - Ai{l3dwgo - d^go) - drgo . (61) 

Now we may consistently require d\gi = 0. In this way we avoid an additional secular 
term and end up with a closed equation for g^. Its solution coincides with (!59|) with the 
first order result 

W(X) = -A (62) 

as expected. Extending the expansion to the next order, we find = as well 
as C2 = 0. The latter is fine because C2 describes the leading deviations from the 
Gaussian shape and in the present case P{W, A) has to remain Gaussian to all orders. 
To understand the former we note that from (l60l) we have for the exact W{X) 

d^W{X) = -{l~e-''''^/') . (63) 
A'" 

As noted already in [5] the second term is non-perturbative and does not show up in 
perturbation theory at any order. 



8. Conclusion 



We have outlined a multiple time-scale perturbation expansion for the systematic 
determination of the work distribution in isothermal quasistatic processes. The 
expansion parameter is the ratio between the relaxation time of the system and the 
typical time scale of the driving. To leading order in this parameter and for typical 
work values the work distribution is Gaussian as expected from the central limit theorem. 
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Higher order terms of the perturbative expansion yield systematic information about the 
deviations from the Gaussian form. However, irrespective of time scale separation the 
expansion breaks down when applied to work values from the tails of the distribution, 
which is again in accordance with the central limit theorem. When considering quantities 
that are not sensitive to the tail of the distribution [16] the Gaussian approximation 
works, of course, well. Often, however, the tails are crucial for averages relevant in 
stochastic thermodynamics [U H] and care must be exercised in using the Gaussian 
approximation. We finally note the pivotal role of the coupling to a heat bath for 
the validity of the Gaussian approximation: in non- isothermal quasi-static processes 
the work distribution may well deviate from a Gaussian form [171 UHl [19] • In the 
adiabatic expansion of an ideal gas, e.g., there are persistent correlations between 
different increments of the work which render the application of the central limit theorem 
impossible [2n|[T9]. 
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